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FOREWORD 


This is a progress report on the research project entitled "Analysis 
and Computation of Internal Flow Field in a Scramjet Engine," for the 
period ended June 30, 1987. Specific attention has been directed to 
investigate "Transient Radiative Energy Transfer in Incompressible Laminar 
Flows." This work was supported by the NASA Langley Research Center 
(Computational Methods Branch of the High-Speed Aerodynamics Division) 
under research grant NAG-1-423. The grant was monitored by Drs. A. Kumar 
and J. P. Drummond of the High-Speed Aerodynamics Division. 



TRANSIENT RADIATIVE ENERGY TRANSFER 
IN INCOMPRESSIBLE LAMINAR FLOWS 


S. N. Tiwari 1 and D. J. Singh 2 
Old Dominion University, Norfolk, Virginia 23529-0247 

ABSTRACT 

Analysis and numerical procedures are presented to investigate the 

transient radiative interactions of nongray absorbing-emitting species in 
laminar fully-developed flows between two parallel plates. The particular 
species considered are OH, CO, C0 2 , and H 2 0 and different mixtures of these 
species. Transient and steady-state results are obtained for the temperature 
distribution and bulk temperature for different plate spacings, wall 

temperatures, and pressures. Results, in general, indicate that the rate of 

radiative heating can be quite high during earlier times. This information is 

useful in designing thermal protection systems for transient operations. 


* Eminent Professor, Dept, of Mechanical Engineering and Mechanics. 
Graduate Research Assistant, Dept, of Mechanical Engineering and Mechanics. 
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NOMENCLATURE 
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t 

T 

T 1 

T fc 

u 

u c 

y 

B 
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band absorptance = A(u,e), cm 

band width parameter, cm -1 

„-l 


-1 


,-l 


correlation parameter, atm x - cm' 
specific heat at constant pressure, kJ/kg-k = erg/gm-k 
Planck's function, (W-cm”^)/cm - ^ 

Planck's function evaluated at wave number w 


o 

e l » e 2 

H li» H l 

k 

K 1 

M 

N 

P 

^R 


„-2 


emissive power of surfaces with temperatures Tj and T^, W-cm 
gas property for the large path length limit 
thermal conductivity, erg/cm-sec-K 

gas property for the optically thin limit distance between plates 
large path length parameter, nondimensional 
optically thin parameter, nondimensional 
pressure, atm 

O 

total radiative heat flux, W/cm 

spectral radiation heat flux, (W-cm _<: )/cm 
wall heat flux, W/cm^ 

1 O 

integrated intensity of a wide band, atm” -cm 
time, sec 
temperature, K 

wall temperature, K; Tj = T w 
bulk temperature, K 
nondimensional coordinate = SPy/A 0 
nondimensional path length = SPL/A 0 
transverse coordinate, cm 
line structures parameter 
nondimensional temperature 



dimensionless bulk temperature 
spectral absorptin coefficient, cnT* 
nondimensional coordinate = y/L = u/u 
density, kg/m^ 
nondimensional time 
wave number, cm - * 

wave number at the band center, cm - * 
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INTRODUCTION 

The radiative energy transfer in participating medium has received 
special attention in recent years because of its applications in the areas of 
the remote sensing, earth's radiation budget studies and climate modeling, 
fire and combustion research, entry and reentry phenomena, hypersonic 
propulsion and defense-oriented research. In most studies involving combined 
mass, momentum, and energy transfer, however, the radiative transfer 
formulation has been coupled mainly with the steady processes [1-20] and the 
interaction of radiation in transient processes has received very little 
attention. However, the transient approach appears to be the logical way of 
formulating a problem in general sense for elegant numerical and computational 
solutions. The steady-state solutions can be obtained as limiting solutions 
for large times. 

The limited number of studies available on the transient radiative 
transfer in gaseous systems [21-31] are reviewed critically in [32-34]. The 
literature survey reveals that the transient behavior of a physical system can 
be influenced significantly in the presence of radiation. 

The goal of this research is to include the nongray radiative formulation 
in the general governing equations and provide the step-by-step analysis and 
solution procedure for several realistic problems. The basic formulations are 
presented in [32] and the specific case of transient radiative exchange in 
nongray gases between two parallel plates is investigated in [33, 34]. The 
objective of this study is to investigate the interaction of transient radia- 
tion in fully-developed laminar flows between two parallel plates. Thus, 
special attention is directed to include the nonsteady phenomenon only in the 
energy equation. Homogeneous as well as nonhomogeneous mixtures of absorbing- 
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emitting species are considered. In subsequent studies, the present analysis 
and numerical techniques will be used to investigate the unsteady flow of 
compressible and chemically reacting species in one- and multi-dimensional 
systems . 
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BASIC FORMULATION 


The physical system considered is the energy transfer in laminar, 
incompressible, constant properties, fully-developed flow of absorbing- 
emitting gases between parallel plates (Fig. 1). The condition of uniform 
surface heat flux for each plate is assumed such that the temperature of the 
plates varies in the axial direction. Extensive treatment of this problem is 
available in the literature under steady state conditions [8, 11]. The 
primary motivation of this study is to investigate the extent of transient 
radiative interaction for high temperature flow conditions. 


For the physical problem considered, the energy equation can be expressed 
as [1] 


» c p $ * 


9T + v il) 

3y j 


3x 


3 2 T 

7 




r 9U \ 2 


divq c 


( 1 ) 


where u and v denote x and y components of velocity, respectively. In 
deriving Eq. (1) it has been assumed that the net conduction heat transfer in 
the x direction is negligible compared with the net conduction in the y direc- 
tion. This represents the physical condition of a large value of the Peclet 
number. By an analogous reasoning, the radiative heat transfer in the x 
direction can be neglected in comparison to that transferred in the y 
direction. If, in addition, it is assumed that the Eckert number of the flow 
is small, then Eq. (1) reduces to 


3 2 T 


1 


aq 


R 


3T . „ 3T . u 3T _ 

IT u 9x 3y ‘ 

where a = (k/pCp) represents the thermal diffusivity of the fluid. 


( 2 ) 


For a steady fully-developed flow, v = 0, and u is given by the well- 


known parabolic profile as 



4 



5 


u = 6 u m (g - 5 2 ); 5 = y/L (3) 

where u m represents the mean fluid velocity. Also, for the flow of a perfect 
gas with uniform heat flux, 3T/3X is constant and is given by 


3T/3X = (2aq w )/(u m L/k) 


(4) 


Now, by combining Eqs. (2) - (4), the energy equation is expressed in 

nondimensional form as 


I 2 - + 12 (5 - C 2 > - 
9t 3{ 


1_ _^R 

q w 


(5) 


where 

T = at/L 2 ; e = (T - T 1 )/(q w L/k) 


By assuming that the initial temperature distribution in the gas is some 
uniform value T 0 = Tj, the initial and boundary conditions for this problem 
can be expressed as 


0 ( 5 , 0 ) = 0 ( 6 a) 

0(0, x) = 0 (1, t) = 0 (6b) 


9 U = 1 / 2 ) = 0 ; 0^(5 = 0 ) = - 0 ? (c = 1 ) ( 6 c) 

It should be noted that all the boundary conditions given in Eqs. ( 6 ) are not 
independent and any two convenient conditions can be used to obtain solutions. 
Also, the initial temperature distribution can be any specified or calculated 
value of 0 ( 5 , 0 ) = f(?) . 

As discussed in Refs. 32-34, the radiative transfer equations are 
formulated for one-dimensional planar systems for many engineering and 
astrophysical applications. For diffuse nonreflecting boundaries and in the 
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absence of scattering, the expression for the total radiative flux is given, 
for an n-band gaseous system, by [1, 8 , 32] 


q R (y) - - e 2 


+ I . 2 , / <J F i„.< z > V ex e [ - 1 <^ z)]dz 

1=1 0 1 1 1 


- / F, (z) k exp [- i < (z-y)]dz}du>. 

y 2 “i W i 2 “i 1 


(7) 


where 


F, (z) = e (z) - e, ; F 0 (z) = e (z) - e 7 

lw. ' lo>.j Zuj oi i 2 a). 


Equation (7) is in proper form for obtaining the nongray solutions of 
molecular species. In fact, this is an ideal equation for the line-by-line 
and narrow-band model formulations. However, in order to be able to use the 
wide band models and correlations, Eq. (7) is transformed in terms of the 
correlation quantities as [1, 7-13, 32] 

( 5 ) - ej - e 2 


♦4 fl oi J i [ 4 u oi (s-s'Hds’ 


/ f 2 oi u') *i ('■-tm') 


( 8 ) 


5 = u/u Q = y/L; 5 ' = u7u Q = z/L; A = A/A q ; 
u = (S/A Q ) py; u Q = (S/A Q ) PL; PS = / d w 


where 
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It should be noted that Fj and in Eq. ( 8 ) represent the values at the 

center of the 1th band and A’(u) denotes the derivative of A(u) with 
respect to u. Upon perform ng the integration by parts, Eq. ( 8 ) can be 
expressed in an alternate form as [32] 


1 r<«> = e l - e 2 

n 
E 

1*1 


+ Z A oi {/* [de u ( 5 ’)/d 5 '] A i [| u Q . ( 5 - 5 ' )]d?' 


+ / [de (5')/d5'] A 1 [f u Qi (s'-ODd?'} 
5 1 


(9) 


A direct differentiation of Eq. (9) provides the expression for the divergence 
of radiatve flux as 


dq R U) 
d? " 


7 A 01 </ 0 [de „ i ( *' )/d «' :l x 


A; [f u 01 .(c-c')]Oc' - / Cde^tc'J/de'] X 

a; [f u 0 ,(5'-c)]dc'} ( 10 ) 

Equations ( 8 ) through (10) are the most convenient equations to use when 
employing the band-model correlations in radiative transfer analyses. 

For the present physical problem, e^ = e£ and F^ = ^ 2 ^-’ Thus, for 
the case of linearized radiation, a combination of Eqs. (5) and ( 8 ) results in 
[ 8 , 32] 

0 . 0 . 3 Ne - 12 ( ? -c 2 ) = j (L/k) Z H u 2 {/* eU', t) x 

55 t 4 i=i n 01 0 

3 1 

Al' [| u oi (5-C')]dc' + / 0(5', t) x 
A]' [f u of ( 5 ’-c)]d ? ') 


(11a) 
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where 

N = (PL 2 /k) K. = (PL 2 /k) e S . ( T) (de /dT) T 
i i=1 i l 

H li = V T > ( de u i /dT > T l * H 1 = 1 = 1 H li 

The dimensionless gas property N characterizes the relative importance of 
radiation versus conduction within the gas under optically thin conditions. 
Also, by combining Eqs. (5) and (10) another form of the transient energy 
equation is obtained as 

e rr -6-12 U-5 2 ) = 1 (L/k) E H . u {/* (30/ac 1 ) X 
55 t ^ i=i u 01 o 

a; [f u oi (s-s')]d 5 ' - jj (89/3?‘) X 

A! [| U oi (c’-5)]d?'} (lib) 

Note again that Ea. (lib) can be obtained directly by integrating the right- 
hand side of Eq (11a) by parts. Quite often, Eq. (lib) is the convenient form 
to use in radiative transfer analyses. 

For flow problems, the quantity of primary interest is the bulk 
temperature of the gas, which may be expressed as [11] 

1 ? 

e b = (T b _T l )/(q w L/k) = 6 / 0 0 (5 >t)U-| )d C (12) 

The heat transfer q w is given by the expression, q w = h c (Tj-T^), where h c is 
the convective heat transfer coefficient (W/cm 2 -K). In general, the heat 
transfer results are expressed in terms of the Nusselt number Nu = h c D^/k. 
Here, represents the hydraulic diameter, and for the parallel plate 
geometry it equals twice the plate separation, i.e., = 2L. Upon 
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eliminating the convective heat transfer coefficient h c from the expressions 
for q w and Nu, a relation between the Nusselt number and the bulk temperature 
is obtained as 

Nu = 2 Lq w /k(T r T b ) = -2/e b (13) 

The heat transfer results, therefore, can be expressed either in terms of Nu 
or 9 b . 

LIMITING CASES AND SOLUTIONS 

Before discussing the method of solution for the general case, it is 
advisable to explore the various limiting cases. Quite often, closed form 
solutions can be obtained for some of these cases. Specifically four limiting 
cases are considered here and these are the steady laminar flows, the case of 
negligible radiation, the optically thin limit, and the large path length 
1 imit. 

Steady Laminar Flow 

For steady-state conditions 30/3t = 0 and Eqs. (11) provide two forms 
of the energy equation for this case. Another convenient form is obtained by 
letting 30/3 t = 0 in Eq. (5) such that 

0 " - 12 U-C 2 ) = (l/q w ) dq R /d? (14) 

By integrating Eq. (14) once and using the conditions that at 5 = 1/2, 
q R (?) and (de/d?) are equal to zero, one obtains 

e' - 2(3? 2 - 2? 3 ) + 1 = q R U)/q w 

A combination of Eqs. ( 8 ) and (15) results in 

0 o n 5 

0 ' - 2(3? - 2? 3 ) + 1=4 (L/k) E H. . u . {/ e(c') x 

L i=l 11 01 0 


(15) 



10 


a; [§ U 0l («-5')]d5' - /* 8(5') x 

a; [| u o ) ( C , -« ) W5') (16) 

It should be pointed out that by combining Eqs. (9) and (15) another useful 
form of the energy equation can be obtained for the steady case. 

For the case of negligible radiation, Eq. (16) reduces to a very simple 
form and utilizing the boundary condition 0(0) = 0, the solution of the 
resulting equation is found to be 

9(5) = 5(2 5 2 - 5 3 - 1) ( 17a ) 

The result for the bulk temperature is found by combining Eqs. (12) and (17a) 
as 

- e b = 17/70 (17b) 

The results provided by Eqs. (17) are useful in determining the extent of 
radiatve contributions. 

Negligible Radiation 

For the case of negligible radiation, N = 0 and both forms of Eq. (11) 
reduce to 

0 55 - 0 T - 1? (5-£ 2 ) (18) 

By employing the product solution procedure, the solution of Eq. (18) can be 
obtained and the result can be expressed in terms of the bulk temperature 
through use of Eq. (12). Alternately, the solution of Eq. (18) is assumed to 
be of the form 


0(5.t) = g(c) + h(g, T ) 


(19) 
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From Eqs. (18) and (19), there is obtained two separate equations as 

g" = 12 U- 5 2 ) (20) 

h - h T = 0 (21) 

Solutions of Eqs. (20) and (21) are found separately from which the 
complete solution for the temperature distribution is obtained as [32] 

0 ( 5 , t) = 5(2 5 2 - C 3 - 1 ) + E C sin (a 5 ) exp(-a t); a 2 = nii ( 22 a) 

n-1 n 

where 

C n = (4/a 5 ) [ ( 12 - 12a 2 + a 4 ) cos(a) - 24], n = 1, 2,... 

Finally, the expression for the bulk temperature is found by using Eq. (12) as 


00 . _ 

9 . = -17/70 + 6 e C [(1/a) + (4/a 3 )] exp(-aS) 


(22b) 


n=l 


Equations (22) are useful in determining the extent of radiative contributions 
for the transient case. 

Optically Thin Limit 

In the optically thin limit, the expression for the bulk temperature for 
the steady case is found to be [11, 32] 


6 b = [1/ ( 3N ) 3 ] {576(3N)" 1/2 (NEXP) - 21.6N 2 + 72N - 288} 


(23) 


where 


NEXP = {1-exp [-( 3N 1/2 ]}/{l + exp[-(3N) 1/2 ]} 

Both forms of the transient energy equation, Eqs. (11a) and (lib), reduce 
to a simplified form in the optically thin limit as [32] 
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Assuming a solution 


h 


55 


e - e - 3Ne = l2(c-c 2 ) ( 24 ) 

55 t 

of the form given by Eq. (19), Eq. (24) is written as 

- h - 3Nh = - g + 3Ng + 12(5 - 5 2 ) (25) 

t 55 


Consequently, 

g" - 3Ng = 12(5 - 5 2 ) (26) 

and h ?? - h T - 3Nh = 0 (27) 

From the solution of Eqs. (26) and (27), the solution for Eq. (24) is obtained 
as 


where 


0 ( 5 , t) = (16/3N 2 ) [sinh (-✓W2)/s1nh(/3R)] 

cosh IM (5-1/2)] + (4/N)(5 2 - 5 + 2/3N) 

+ z C sin (a 5 ) exp [- ( 3N + a 2 ) t]; a = nir (28a) 

n=l n 


C = 0 , for n even (28b) 

n 

= 32 ( 3N + a 2 )/(3N 2 a 3 ) 

+ 2a/[3N 2 (3N + a 2 )] , for n odd (28c) 

By combining Eqs. (12) and (28), the expression for the bulk temperature is 
obtained as 

9 b = 6 { ( 1 6/3N 2 ) [si nh ( -/3N/ 2 ) /si nh ( /3N ) ] x 
[ ( 1/3N) cosh(/3N/2) - (4 + /3N) (3N)‘ 3/2 x 
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sinh(/3R/2)] + (4/N) [-1/30 + 1/(9N)] 

+ z C (1/a + 4/a^) exp[-(3N + a 2 )x]} (29) 

n=l n 

Large Path Length Limit 

In this limit, the steady-state energy equation, Eq. (16), reduces to [8, 
11, 32] 

o o 1 

e' - 2(3r - 2r) + 1 = M / 0(5 1 ) de 7 ( c - %') (3°) 

where 

n 

M = Hj L/k = (L/k) i A Qi (de ui /dT) T ^ 

The nondimensional parameter M constitutes the radiation-conduction 
interaction parameter for the large path length limit. Equation (30) does not 
appear to possess a closed form solution; a numerical solution, however, can 
be obtained easily. 

In the large path length limit, the transient energy equations, Eqs. 
(11a) and (lib) , reduce to 

0 — 0 - 3Ne - 12(5 - 5 2 ) = - M [J 5 0 ( 5 ', t) dc'/U - 5*) 2 

55 T 0 

+ / 0(5' ,T)d5'/(?'-5) 2 ] ( 31a ) 

5 

0-0 - 12 (5 - 5 2 ) = M / (30/35' )d5'/U-5' ) (31b) 

55 T ■ 0 

Since | (5-5' ) 2 1 = | (5' -5) 2 | , Eq. (31a) can be written as 

5 1 ? 

0 0 . 3Ne - 12(5 - 5 2 ) = - M / 0 ( 5 ', t) ds 7(5 - 5') 

£5 t n 


(31c) 
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Through integration by parts, Ecu (31c) can be expressed as 

o 1 

e - e - 3Ne - l2(c - c) = M / (ae/ac' )dc7(c-c' ) ( 31d ) 

55 T 0 

Equations (31a) - (31d) represent different forms of the governing equations 
in the large path length limit. With the exception of the term (-3Ne) on 
the left-hand side, Eq. (31d) is identical to Eq. (31b). Since N represents 
the radiation-conduction interaction parameter only in the optically thin 
limit [ 8 ], it should not appear in the governing equation for the large path 
length limit. Thus, Eq. (31b) is the correct equation to use for solution in 
the large path length limit; the solution of this equation is obtained by 
numerical techniques. 


METHOD OF SOLUTION 

The solution procedures for both steady and unsteady cases are presented 
in this section. In principle, the same numerical procedure applies to both 
the general and large path length limit cases. 

Steady-State Solutions 

The general solution of Eq. (16) is obtained numerically by employing the 


method of variation of parameters. For this, a polynomial form for eU) is 
assumed in powers of 5 as 

6(c) = J a m ? m (32) 

m =0 

By considering a five term series solution (a quartic solution in 5 ) and 
satisfying the boundary conditions 0 ( 0 ) = 0 1 ( 1 / 2 ) = 0 and e'( 0 ) = -e'(l), 

one obtains 

0(C) = ajU - 2c 3 + C 4 ) + a 2 (?2 ' 2 ^ + ^ (33) 
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e'U) = a ! ( 1 - 6(- 2 + H 3 ) + a 2 (2? ' + 4?3 ) (34) 

A substitution of Eq. (34) in Eq. (16) results in 

«j(l - 6 e 2 + 4? 3 ) ♦ a 2 (2 5 - 6{ 2 ♦ 4| 3 ) 

- 2(3 { 2 - 2{ 3 ) *1=| (L/k) l H X1 u 0 , x 

(/' 6(5') A; [|u o1 ( 5 - 5')ld 5 ' 

- / e(e') A! [| u ,({' - 5)] dc'J (35) 

5 

where expressions for 0 ( 5 ') are obtained from Eq. (33). 

The two unknown constants aj and a 2 in Eq. (35) are evaluated by 

satisfying the integral equation at two convenient locations U=0 and 5 = 1/4 
in the present case). The entire procedure for obtaining aj and a 2 is 
described in [32]. With known values of and a 2 , Eq. (33) provides the 
general solution for 0 ( 5 ). The expression for the bulk temperature is 

obtained by combining Eqs. (12) and (33) as 

e b = (1/70) ( 17a j + 3a 2 ) (36) 

It should be noted that for the case of no radiative interaction a 2 = 0 and a± 

= -1, and Eq. (36) gives the result of Eq. (17). 

The governing equation for the large path length limit is Eq. (30). For 
this equation also the solution is given by Eqs. (33) and (36) but the values 
of a's are completely different in this case [32]. 
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Transient Solutions 

The governing energy equations for the transient case are Eqs. (11a) and 
(lib). General solutions of these equations are obtained also numerically by 
employing the method of variation of parameters. For the present problem, a 
polynomial form for 6 ( 5 , t) is assumed as 

e(e.t) = z a (t) < 37 ) 

m=0 

For a quadratic temperature distribution in ? (with time dependent 
coefficients), Eq. (37) is written as 

eU,x) = 3 q(t) + aj(t) 5 + a 2 ( T ) ^ (38a) 

By using the boundary conditions 0(0 , t) = 0 and 0^(5 = l/2) = 0, this reduces to 

0 ( 5 , t) = g(x) ( 5 -c 2 ) (38b) 

where g(t) represents the time dependent coefficient. Consequently, 

0 (C, t) = g(x) (1-25); 9 5? (5,t) = -2g(x); 0 t (?,t) = (C5 2 ) 9'(t) (39) 

Also, a combination of Eq. ( 6 a) and (38b) yields the initial condition 

e(c,o) = g(0) = 0 (40) 

Note that essential boundary conditions are used already in obtaining the 
solution represented by Eq. (38b). 

By employing Eqs. (38b) and (39), Eqs. (11a) and (lib) are transformed in 
alternate forms which are expressed in a compact form as 

J^?), 

J 2 U) 


g' ( t) + 


g( t) + 12 = 0 


(41) 
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where JjU) and J 2 (g) are defined in [32]. The function ^( 5 ) is used for 
solution of Eq. (11a) and J 2 U) is used for solution of Eq. (lib). The 
solution of Eq. (41) satisfying the initial conditions of Eq. (40) is given by 

g( t) = ^J||y {exp C-J(c)t]-1} (42) 

The temperature distribution given by Eq. (38b) can be expressed now as 

9 ( 5 , t) = -j||y {exp[-J(?)x]-l}(c-5 2 ) (43) 

The expression for the bulk temperature is obtained through use of Eq. (12) as 

9 , = 72 / [U-S 2 ) 2 /JU)]{exp[-J( 5 )x]-l} (44) 

D 0 

Note that in Eqs. (42)-(44), J(§) becomes JjU) for solution of Eq. (11a) 
and J 2 (5) for solution of Eq. (lib). 

For a quartic solution in 5 , Eq. (37) gives the result similar to Eq. 
( 33 ) which for the transient case is expressed as 

9(5, t) = g(x) U - 25 3 + € 4 ) + h(x) (C 2 - 2 ? 3 + 5 4 ) (45) 

By substituting Eq. (45) into Eq. (11a), one obtains 

x g'(x) + J 3 U) 9 (t) + yh'(x) + J 4 U) h ( t) = - z (46) 

where 

x = (5 - 2c 3 + 5 4 ); y = U 2 - 2? 3 + c 4 ); 
z = 12( 5 - c 2 ) 

and functions J^U) and J 4 U) are defined in [32], Equation (46) 

constitutes one equation in two unknowns, namely g(x) and h(x). However, 
since the equation is linear in x, the principle of superposition can be used 
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to split the solution into two equations as 

x g ’ ( t) + J 3 ( 5 ) g(x) = -z/2 ( 47 ) 

yh'(x) + J 4 (c) h(x) = -z/2 (48) 

The initial condition for this case can be written as 

9 ( 5 , 0 ) = g( 0 ) U - 25 + § 3 ) 

+ h ( 0 ) (c 2 - 2c 3 + c 4 ) - 0 (49a) 

Consequently, 

g ( 0) = 0; h(0) = 0 (49b) 

The solution of Eqs. (47) and (48) satisfying the appropriate initial 
condition of Eq. (49b) is given respectively as 

g(i) = Cz(c)/ 2 J 3 (c)] {exp[-J 3 U)T/xU)]-l} (50) 

h ( t ) = [z(?)/ 2J 4 U)] {exp[-J 4 ( 5 )x/yU)]-l> ( 51 ) 

By substituting Eqs. (50) and (51) into Eq. (45), the expression for the 
temperature distribution is obtained as 

9(5, x) « [6(5-5 2 )(5-25 3 +5 4 )/0 3 (c)] x 

{exp[-J 3 ( 5 )x/x(c)]-l} + [6(c-c 2 ) x 

(5 2 -25 3+ 5 4 )/J 4 (5)] {exp[-J 4 (5)t/y(5)]-l> 

The bulk temperature in this case is given by 

0 K = 36 / [(5-5 2 )(5-25 3 +5 4 )/0 3 (?)] x 
° 0 


( 52 ) 
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{expC-J^(5)T/x(c)]-l}d 5 + 36 / [( 5 -C 2 ) x 

J 0 

(C 2 -2? 3+ ? 4 )/J 4 (c)]{exp[-J 4 (§)T/y(5)]-l}d?. (53) 

where x and y are defined in Eq. (46). 

By substituting Eq. (45) into Eq. (lib), there is obtained 

xg' + J 5 (?) g(t) + yh + J 6 U) h( T ) = -z (54) 

where again x,y,z are defined in Eq. (46) and functions Jg( 5 ) ar| d ^g(s) are 
defined in [32]. The solution procedure for this equation is identical to 
that for Eq. (46) and the results for temperature distribution and bulk 
temperature are given respectively by Eqs. (52) and (53) with J 3 replaced by 
J 5 and J 4 by Jg. 

In the large path length limit, the two applicable governing equations 
are Eqs. (31b) and (31d). The solutions of these equations can be obtained 
from the general solutions by evaluating the integrals in J function in the 
large path length limit [32]. 

PHYSICAL CONDITIONS AND DATA SOURCE 

As discussed in [32-34], four specific absorbing emitting species were 
selected for an extensive study; these are CO, CO 2 , OH and H 2 O. The species CO 
was selected because it contains only one fundamental vibration-rotation (VR) 
band and all spectral information are easily available in the literature. It 
is a very convenient gas to test the numerical procedure without requiring 
excessive computatinal resources. Species OH and H 2 O are the primary 
radiation participating species for the pressure and temperature range 
anticipated in the combustor of the scramjet engine. Species CO 2 , and 
combinations of C0 2 and H 2 0 are important absorbing-emitting species in many 
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other combustion processes. Different mixtures of various species (such as 
C0 2 + H 2 0, OH + H 2 0, and CO 2 + H 2 0 + OH) were selected for parametric 
studies. Thermophysical properties of these species are given in [33] for 
different temperatures. 

In radiative transfer analyses, it is essential to employ a suitable 
model to represent the absorption-emission characteristics of specific species 
under investigation. Several line-by-line (LBL), narrow-band, and wide-band 
models are available to model the absorption of a VR band [7-12]. However, it 
is often desirable to use a simple correlation to represent the total 
absorption of a wide band. Several such correlations are available in the 
literature [7-12]. The relative merits of these correlations are discussed in 
[12]. In this study, the correlation proposed by Tien and Lowder [7] is 
employed and this is given by 

S(u) = In {uf ( e ) [— + + (55) 

where 

f ( 6 ) = 2.94[1 - exp(- 2.60 e)] 

and 6 represents the line structure parameter. 

The spectral information and correlation quantities and thermodynamic and 
transport properties needed for different species were obtained from Refs. 7-9 
and 33. Spectral information and correlation quantities for important 
infrared bands are provided in Appendix A. The specific VR bands considered 
for each species are: CO (4.7 y fundamental), OH (2.8 y fundamental), C0 2 (15 y , 

4.3y and 2.7p), and H 2 0 (20p, 6.3p, 2.7p, 1.87 u , and 1.38y). 

In a mixture of several species, spectral lines and bands overlap in 
certain spectral regions. The total absorptance in such regions cannot be 
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calculated simply by adding the contributions of different bands and 
corrections should be made to account for the partial overlapping. If line- 
by-line or narrow band models are employed in the general formulation of the 
physical problem, then there is no need for such corrections [12, 35, 36]. 
The solution of LBL formulation, however, requires considerably large 
computational resources. Use of narrow band models offers some computational 
relief but certain spectral information needed are not available for many 
species for temperatures higher than about 600 K. A relatively easier 
procedure (called the block method) is suggested by Edwards [13] and is useful 
in calculating the total emissivity of a mixture of several species. Another 
method suggested by Penner and Varanasi [37] is probably the most convenient 
method to use in the frame work of the radiative flux formulatin expressed in 
terms of the wide-band model absorptance and correlations. 


For a homogeneous path, the total absorptance of a band is given by 


oo 


A(y) = / [1 - exp(- k y)] du (56) 

o “ 

where both < and have units of cm If in a spectral range Au-, there 
u 1 

are contributions from bands of different species, then for a homogeneous path 
the transmittance is given by 


N 

x. = exp(- e < y) 
Au i . j=l “ij 


(57) 


where N represents the number of participating species in the gaseous mixture. 
Consequently, Eq. (56) can be expressed as 

N 

A i = J [! - ex P(- * V. 

Au>, J=1 1 J 


y)] du> 


(58) 
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If two bands of different (or same) species are occupying approximately the 
same spectral range Au>.j , then Eq. (58) reduces to 

A. = / {1 - expfic + k )y]} don (59) 

1 “11 “12 


By employing the relations for the exponentials, Eqs. (59) can be expressed as 


[37] 



A, * / [1 

A(jj * 

1 

- exp(- k y)] da) + / [1 - exp(- )] dco - aA 

“i 1 Au).j 12 

(60a) 

or 

A = A i + A 2 - aA 

(60b) 

where 




aA - / {[1 - exp(- y)] x 

Aa).j il 


[1 - exp(- y)]} do) (60c) 

Use of Eqs. (60) has been made by Felske and Tien [38] and Tien and Lee [39] 
to calculate absorptances of homogeneous and nonhomogeneous mixtures of CO 2 
and H 2 0 in the 2.7 y region for different pressure and temperature 

conditions. A similar procedure is used in this study to account for the 

overlapping effects of different species. 

For the physical problem considered, the dependent variables are e and 
e b (or Nu) and independent variables are t and 5 . The parameters, in 

general, are T w , P, and L. The large path length and optically thin limits 

are characterized respectively by parameters M and N. Radiative and thermo- 
physical properties of participating species are evaluated at different 

specified pressures and temperatures. 
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RESULTS AND DISCUSSIONS 

Extensive results have been obtained for variation of e and e b for 
different conditions and most of these are available in Appendix B. The 
computer program used for numerical solutions is provided in Appendix C. 
Selected results are presented here to compare solutions of quadratic and 
quartic formulations and demonstrate the variation of e with t and 5 and 
of e b with L for single component systems and homogeneous mixtures. For all 
results presented here, a lower value of e (or e b ) in the figures indicates 
a higher value of temperature in the medium; this, in turn, implies a rela- 
tively higher ability of the gas to transfer radiative energy. 

Results of quadratic and quartic formulations are compared in Figs. 2-5 
for different species. The centerline temperature variations with nondimen- 
sional time are compared in Fig. 2 for P = 1 atm, T w = 500 K, and L = 5 cm. 
The results show that the steady-state conditions are reached at an eariler 
time for H 2 0; and this is followed respectively by C0 2 , CO, and OH. 
Significant differences are noted between the quadratic and quartic solutions 
for larger times. For the conditions of results presented in the figure, the 
difference are found to be greatest for OH and lowest for H 2 0. The results 
for 9 versus 5 are compared for P = 1 atm, L = 10 cm, and T w = 500 K in 
Fig. 3 and T w = 1,000 K in Fig. 4. The results demonstrate that considerable 
differences in solutions can occur at different locations in the channel and 
that the differences are larger for the lower wall temperature. Results for 
e b versus L presented in Fig. 5 show that the differences in two solutions 
are relatively larger at lower plate spacings and that quartic results 
approach the correct limiting solution for the case of no radiative interac- 
tion (9 b = -0.243). From the results presented in Figs. 2-5 and in Refs. 33 
and 34 it is concluded that while quadratic and quartic solutions are 
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Figure 2. Comparison of quadratic and quartic results for the centerline 
temperature variations with time; T = 500 K, P = 1 atm. 
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Figure 3. Comparison of quadratic and quartic results for temperature 
variations across the duct; P = 1 atm, L = 10 cm, and 




Figure 4. Comparison of quadratic and quartic results for temperature 
variations across the duct; P = 1 atm, L = 10 cm, and 
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Figure 5. Comparison of quadratic and quartic results for the bulk 
temperature; P = 1 atm and T w = 500 K. 
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identical for radiative equilibrium and radiation with conduction cases, they 
differ significantly for the case of combined conduction, convection and 
radiation. As such, all other results for this study were obtained by using 
the quartic formulation. 

The results for temperature variations with nondimensional time are 
presented in Figs. 6-9 for various species and for different physical 
conditions. The centerline temperature distribution for general and limiting 
cases are illustrated in Fig. 6 for T w = 500 K, P = 1 atm, and L = 5 cm. The 
results show that the steady-state conditions are reached at about t = 0.5 
for all species. As noted earlier, these results also demonstrate that H 2 0 is 
a highly radiation participating gas as compared to C0 2 , CO, and OH. In 
comparison to other species considered, OH takes relatively longer times to 
reach the steady state and is least effective in transferring the radiative 
energy. For the specified physical conditions, the large path length solu- 
tions are closer to the general solutions and optically thin solutions provide 
higher rate of energy transfer. For all species, optically thin solutions 
reach the steady state faster than other solutions. Also, in the optically 
thin limit, C0 2 is more effective in the radiative transfer process than other 
species. The reasons for such trends are given in [8, 33, 34]. The results 
for T = 500 K and 1,000 K are compared in Fig. 7, and they simply indicate 
that the rate of energy transfer is higher at the higher temperature. The 
results for 5 = 0.25 and 0.5 are compared in Fig. 8 for T w = 1,000 K, P = 1 
atm, and L = 5 cm. It is noted that the rate of energy transfer is higher at 
earlier tines and at locations closer to the wall. The centerline temperature 
variations for L = 5 cm and 10 cm are illustrated in Fig. 9 for T w - 1,000 K 
and P = 1 atm. As would be expected, the rate of energy transfer is seen to 
increase with the increasing path length. 
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Figure 6. Comparison of general and limiting results for the centerline 
temperature variations with time; T w = 500 K, P = 1 atm, and 
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Figure 7. Comparison of centerline temperature variations with time for 
T w = 500 K and 1,000 K; P = 1 atm and L = 5 cm. 




Figure 8. Comparison of temperature variations with time for = 0.25 and 
0.5; T = 1,000 K, P = 1 atm, and L = 5 cm. 




Figure 9. Comparison of centerline temperature variations with time for 
L = 5 cm and 10 cm; T w = 1,000 K and P = 1 atm. 
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Fiqure 11. Comparison of temperature variations across the duct for OH 
T = 500 K, P = 1 atm, and L = 10 cm. 
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Fiqure 14b. Comparison of temperature variations across the duct for H 2 0 


T = 1,000 K, P = 1 atm, and L = 5 cm. 
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The bulk temperature results as a function of the distance between the 
plates are presented in Figs. 15-19 for different times. General as well as 
limiting solutions are illustrated in these figures. As mentioned earlier, 
the limiting value of = - 0.243 corresponds to negligible radiation. For 
all species, the results presented in the figures are for a pressure of one 
atmosphere. However, for any particular L, the large path length results 
essentially represent the limiting solutions for high pressures. The results, 
in general, demonstrate that the rate of energy transfer is higher at earlier 
times, the effect of radiation increases with increasing plate spacing, and 
the radiative transfer is more pronounced at the higher wall temperature. 

General as well as limiting solutions for the bulk temperature are 
illustrated in Figs. 15-18 for individual species. It is seen that for all 
species the general solutions for t = 0.5 and 1.0 are essentially the same 
for all plate spacings, and the large path length results are valid for 
spacings greater than L = 10 cm for all times. The results for OH are 
presented in Figs. 15a and 15b for T w = 500 K and 1,000 K, respectively. It 
is noted that optically thin results provide the correct limiting solutions 
for plate spacings upto L = 3 cm for all times. The results for T w = 500 K 
show only slight difference between general and large path length solutions 
for t = 0.5, and no significant difference was noted at earlier times (Fig. 
15a). This, however, is not the case for the results presented in Fig. 15b at 
T w = 1,000 K. This trend in results for OH was noted also in [33, 34]. The 
results for CO are illustrated in Fig. 16 for T w = 500 K and P = 1 atm. In 
this case, the optically thin solutions are seen to be valid only upto L = 1.5 
cm. For T w = 1,000 K and P = 1 atm, the results presented in Fig. 17 for C0 2 
and in Fig. 18 for H 2 0 show the same general trend but the extent of radiative 
interaction is entirely different. For C0 2 , there is a considerable 
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difference in general and optically thin solutions for all times. For H 2 0, 
however, the optically thin results are closer to the general solutions for 
spacings upto L = 2.5 cm. The results presented in Figs. 15-18 clearly reveal 
that for a fixed spacing between the plates the rate of radiative heating will 
be considerably higher at earlier times than at the steady state. Thus, in a 
particular physical system, the extent of radiative heating can be very 
intense during the initial stages of operation. It is also important to note 
that both the optically thin and large path length results overestimate the 
influence of radiation. Since these solutions can be obtained with less 
numerical complications, they can be utilized to assess whether or not, for a 
given gas, the interaction of radiation is going to be important. 

A comparison of the general band absorptance results for the four gases 
is shown in Figs. 19a and 19b for a pressure of one atmosphere. The results 
clearly demonstrate the relative ability of the four species for radiative 
transfer at different path lengths. For a plate spacing of greater than L = 
3.0 cm, the results show the same trend as noted in Figs. 2-10. For lower 
plate spacings and relatively higher temperatures, however, C0 2 shows a 
significantly higher ability for radiative transfer than other species. This 
is a typical distinguishing feature of the C0 2 under optically thin conditions 
[8, 33, 34]. 

For steady-state conditions, bulk temperature results were obtained for 
mixtures of different absorbing-emitting species under various conditions and 
some of these are presented in Figs. 20-27. The relative amount of each 
species in the mixture and the relative ability of the species for radiative 
transfer in a given physical condition determine the extent of radiative 


interaction. 



42 


r> 



E 

u 




Figure 15a. Variation of bulk temperature with plate spacing for OH; 



X 

o 

tfi 

O 

© 



in 

© 


i 

i 

i 

_L 

O 

© 


I 


I 




-LARGE PATH LENGTH T w *1,000 K 

- OPTICALLY THIN P = 1 atm 


-0 b = O.243 



M bj/ |^x — M j.)= q e 


Figure 15b. Variation of bulk temperature with plate spacing for OH 
T = 1,000 K and P = 1 atm. 
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Figure 16. Variation of bulk temperature with plate spacing for CO; 
T,., = 500 K and P = 1 atm. 




Figure 17. Variation of bulk temperature with plate spacing for CO 
T.. = 1.000 K and P = 1 atm. 
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19a. Comparison of bulk temperature results for T = 1,000 K and 
P = 1 atm. w 
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As mentioned before, in many combustion processes involving fossil fuels, 
the predominant products of combustion are carbon dioxide and water vapor. 
Extensive studies have been conducted in the literature to determine the total 
emissivity of C0 2 and H 2 0 for homogeneous and nonhomogeneous conditions. 
However, only limited studies are available involving mixtures of C0 2 and H 2 0 
for realistic and important physical conditions. The bulk temperature results 
for three different mixtures of C0 2 + H 2 0 are presented in Figs. 20-23 for 
different temperatures and pressures. As would be expected, the results, in 
general, show that the extent of radiative interaction increases with 
increasing temperature, pressure, and path length. It is seen that the 
radiative interaction is stronger for higher amount of H 2 0 in the mixture. 
However, for the optically thin conditions, the radiative contribution is seen 
to increase with increasing amount of C0 2 at relatively higher temperatures. 

The bulk temperatures results for mixtures of OH and H 2 0 are illustrated 
in Figs. 24-26 and for mixtures of C0 2 , H 2 0, and OH in Fig. 27 for different 
pressures and temperatures. All these results clearly demonstrate that the 
radiative ability of a gaseous mixture essentially depends on the amount of 
highly radiation participating species in the mixture. For example, a 
comparison of results presented in Figs. 24 and 27 for P = 1 atm indicates 
that the rate of radiative transfer is significantly higher with the inclusion 
of 20% C0 2 in the mixture of OH and H 2 0. It is further noted that OH becomes 
a highly radiation participating species at higher temperatures and pressures. 



0.25 C0 2 ♦ 0.75 H 2 0 P = 0.1 atm 
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Figure 20. Bulk temperature results for CO2 + ^0 for P = 0.1 atm 
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Figure 23. Bulk temperature results for COg + H 2 O for T w - 2,000 K. 
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Figure 26. Bulk temperature results for OH and H 2 0 for T w = 2,000 K. 
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CONCLUDING REMARKS 

Analytical formulations and numerical procedures have been developed to 
investigate the transient radiative interaction of absorbing-emitting species 
in laminar fully-developed flows. Extensive results have been obtained for 
OH, CO, C0 2 and H 2 0 for different physical conditions. Illustrative results 
for the temperature distribution and bulk temperature are presented for 
different pressures and wall temperatures. In these results, a lower value of 
temperature implies a relatively higher ability of the gas to transfer 
radiative energy. 

Comparative results of quadratic and quartic formulations confirm the 
need to use the quartic formulation in the numerical procedures for the case 
of combined conduction, convection, and radiation. The results, in general, 
demonstrate that the steady-state conditions are reached at about t = 0.5 
for all species, H 2 0 is a highly radiation participating species (as compared 
to C0 2 , CO, and OH), the rate of energy transfer is higher at earlier times 
and at locations closer to the wall, differences between the limiting and 
general solutions are small at earlier times, the effect of radiation 
increases with increasing plate spacing, and the radiative transfer is more 
prounounced at higher wall temperature and pressure. Similar conclusions can 
be drawn from the results presented for various mixtures. The results clearly 
show that for a given physical condition the radiative ability of a gaseous 
mixture depends essentially on the amount of highly radiating species in the 
mixture. From the results presented in this study, the extent of total 
heating can be determined for different times. This information is essential 
in designing thermal protection systems for operations during the initial 
stages of intense heating. 
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APPENDIX A 


SPECTRAL INFORMATION AND CORRELATION QUANTITIES 
FOR IMPORTANT INFRARED BANDS 


Spectral information and correlation quantities for important infrared 
bands are available in [7-9]. The exponential band model correlation 
quantities for important bands of CO, C0 2 , H 2 0, and CH 4 are given in Table 
A.l. For some species, revised data are now available in [13]. 

In developing correlation quantities for different bands, it is 
important to have information on the integrated band intensity S(T q ) and the 
band width parameter A (T). The band intensities of many molecular bands 
are available in the literature [7-9, 13, 33] but this is not the case with 
the band width parameter. Different relations for A q (T) are proposed by 
Edwards et al. in references cited in [9], From a critical evaluation of 
different relations available in the literature for A Q (T) and after personal 
communications with Dr. D.K. Edwards, it was decided to adopt the following 

two relations: 



2 (3j (4kTBjl/2 
4 1 he 


(A.l) 


A (T) = 0.9 r 2 (1) ( 
0 4 


2kTB-j 1/2 
he 


(A. 2) 


where B is the rotational constant of lower level, c is the speed of light, 
h is the Planck's constant and k is the Boltzmann's constant. The value of 
B is different for different molecules. Equations (A.l) and (A. 2) are 



h-C 


essentially the same except Eq. (A.l) will result in a coefficient of 0.707 

instead of 0.9. It is suggested to use Eq. (A. 2) in determining approximate 

relations for A (T) for all molecules whose values are not available in the 
o 

literature. Information on the rotational constants for different 
molecules is available in the literature (for example, see Ref. 33). It is 
suggested to use the equivalent value for the rotational constant of 
polyatomic molecules. 

It should be noted that in Eqs. (A.l) and (A. 2), B, c, h, and k are 
constants and do not depend on the temperature. Thus, Eq. (A. 2) may be 
expressed as: 


A Q (T) = CONST (T) 1/2 (A. 3a) 

where 

CONST = 0.9 T 2 (2j f 2kBe 1 1/2 (A. 3b) 

4 he 

and B g represents the equivalent rotational constant. By evaluating Eq. 

(A. 3a) at a reference temperature T the value of A Q (T ref ) can be 
determined and, therefore, Eq. (A. 3a) may be expressed alternately as: 

V T > * V T ref> < T/T ref> 1/2 . (A ‘ 4) 

Equation (A. 4) is a convenient form to compare its results with experimental 
values . 

By noting that r 2 (2) = [r (3/4) ] 2 and substituting values for c, h, 

4 

and k, Eq. (A. 3) can be expressed as: 

A (T) = 1.59313 (B e T) 1/2 


(A. 5) 
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where A Q and B g have units of cm 1 and T is in degrees Kelvin. For a 

particular gas, A (T ,) can be obtained from Eq. (A. 5) and then Eg. fA.4; 

K o ref 

can be used to determine A q (T) at other temperatures. For example, for CO 

the rotational constant is 1.931 cm' 1 and at a reference temperature of 

300 K, A (T , = 300 K) = 38.344 cm" 1 . This compares very well with the 
o ref 

experimental value of 38.1 given if Refs. 7 and 9 and presented in Table 
A.l. Similarly, for the 4.3 u band of C0 2> the equivalent rotational 
constant is 0.3906 cm' 1 and, therefore, A Q (T ref = 300 K) = 17.246 cm" 1 ; the 
experimental value of 19.9 cm' 1 given in Table A.l for this band is slightly 
higher. 

Spectral Information for OH 

For the fundamental band of OH, the following information is obtained 
from Ref. 33 : 


Band center, w. = 3570 cm ^ 

-2 -1 

Band strength at STP, S(T q ) = 110 cm atm 
Also, from Ref. 33 the information on equilibrium rotational constant for OH 


A 2 I + ♦ 17.355 cm' 1 
X 2 n i -*• 18.871 cm' 


is obtained as: 
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Thus, it is suggested to use a value for tne equilibrium rotational constant 
for OH as B = 18 cm Using this, a value of A q ( T r efr = 30U K) is found 

as 

A q (T ref = 300) = 1.59313 (B e T ref ) 1/2 = 117.0707. (A. 6) 

Considering the value of A Q (T ref ) = 117 cm -1 , the relation for A q (T) for OH 
is given as: 

1 /? 

A 0 (T) = 117 (T/300) . (A. 7) 

By knowing w , S (T Q ) , and A o (T), other required spectral information for 
OH can be evaluated. 



Table A.l Exponential Band Model Correlation (fciantities 



♦♦Notes on Table A.l. (See next page) 
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]. Correlation quantities are based on the results of reference [9]. The 
intensity of the band marked with * was taken from the reference [7], 

2. Notations: K^T) « (T/300) 172 , K 2 (T) = (300/T) 372 , 

2 (T)/K 1 (T)3 x lo" 3 , fi 2 « 2 (T)/[4> 2 (T)K 1 (T)] 

h = 6.625 x 10 ~ 27 erg-sec, C = 2.998 x 10^ cm/sec, 
k = 1.380 x 10 ^ erg/K , hC/k = 1.44 cm - K 

3. Temperature range: 300 K < T < T max . For CO, T max = 1800 K. 

For C0 2 , T max = 1400 K. For H 2 0, T max = 1100 K. For CH 4 » 


^max 


= 830 K. 


4. For CO, oj = 2143 cm * and 

♦jtT) = [15.15 + 0.22 (T/T q ) 3/2 ] [1-exp (-hCw/kT)], T Q = 100 K 

5. For C0 2 , = 1351 cm *, u> 2 = 667 cm Wj = 2396 cm 


4> 2 ( T) = (1 - exp [( - hC/KT ) (“ 1 + 

{[1 - exp ( - hCoi^/kT)] [1 - exp (-hCio^/kT) ]} 1 , (T ) = 1 + 
3/2 


0.053 (T/100) 



6. For H 2 0, Wj = 3652 cm" 1 , <o 2 = 1595 cm" 1 , « 3 = 3756 cm -1 

* v l v 2 v 3 (T) = U - exp [ - hC(v 1 Wj + v 2 u> 2 + v 3 « 3 )/kT]} x 

{[1-exp ( -hCwj/kT)] [l-exp(-hCo) 2 /kT] [l-exp( -hCu) 3 /kT] > 

- 1/2 

(T) = exp [ - 17.6 (T/100) ] 



APPENDIX B 


RESULTS OF THE QUARTIC FORMULATION 

Extensive results of the quartic formulation are presented in this 

appendix. The order of presentation of results is e vs. t, e vs. c, and 

0 . vs. L. The order of results for the four species considered is OH, CO, 

D 

CO^, and H^O. Information on other physical conditions is provided on each 

figure. The notations used in these figures correspond to the notations 

used in the main text according to the following relations: 

LARGE PATH = LARGE PATH LENGTH 

OPT THIN = OPTICALLY THIN 

TBULK = 6 b = (T w - T b )/(q w L/K) 

THETA = 6 = (T - T)/q L/K) 
w w 

T = T 
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Figure B.l. (continued) 
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Figure B.l. (continued) 
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Figure B.l. (continued) 
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Figure B.2. Results for temperature variation with time (0 vs. t) for CO 
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Figure B.2. (continued) 
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Figure B.2 (continued) 
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Figure B.2 (continued) 
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Figure B.2 (continued) 
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Figure B.2 (continued) 
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The results for temperature variations within the plates are presented in 
Figs. 10-14 for different conditions. Since the temperature profiles are 
symmetric, most of these results are illustrated only for C = 0 to 
S = 0.5. The stady-state solutions for various species are compared in Fig. 
10 for P = 1 atm, L ■ 10 cm, and different wall temperatures. Results again 
demonstrate the relative importance of different species for energy transfer. 
It is noted that the results for CO 2 at T^ = 500 K are exactly the same as for 
OH at T w = 1,000 K, and the results for CO and H 2 O are about the same at T w = 
1,000 K and 2,000 K. This obviously is a coincidence for the physical case 
considered. The general and limiting solutions for OH, CO, C0 2 , and H 2 0 are 
presented, respectively, in Figs. 11 through 14a for T w = 500 K, P - 1 atm, 
and L = 10 cm. These results also show that the steady-state conditions reach 
earlier for H 2 0 and C0 2 than for CO and OH. In general, the differences 
between limiting and general solutions are found to be small at earlier times; 
the maximum difference occurs at the steady-state conditions. In each case, 
the large path length solutions are closer to the general solutions, but the 
optically thin solutions are found to deviate considerably. This is because 
for P = 1 atm and L = 10 cm, the pressure path length is sufficiently high and 
the optically thin limit is not the correct limit for the physical case 
considered. It should be noted that for CO, CO 2 , and H 2 O, the optically thin 
solutions are identical for t > 0.05. Additional results given in Appendix B 
reveal that the differences between general and large path length solutions 
are insignificant at higher temperatures for all species. The general 
solutions for the temperature variation across the entire duct are illustrated 
in Fig. 14b for H 2 0 for the physical conditions of T w = 1,000 K, P = 1 atm, 
and L = 5 cm. The parabolic nature of the transient profiles is clearly 
evident and, in this case, the steady-state is reached at t > 0.5. 
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Figure B.3 (continued) 
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Figure B.4. Results for temperature variation with time (e vs. t) for H 2 0. 
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re B.4 (continued) 
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Figure B.4 (continued) 
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Figure B.4 (continued) 
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Figure B.6. Results for temperature vatiation across the duct (e vs. ?) for CO 
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Figure B.6 (continued) 
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Figure B.7. Results for temperature variation across the duct (e vs. c) for CO 
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Figure B.7 (continued) 
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Figure B.8 (continued) 
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Figure B.9. Results for bulk temperature variation with plate spacing (ej, vs. L) for OH 
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Figure B.9 (continued) 
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Figure B.9 (continued) 
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Figure B.9 (continued) 
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Figure B.9 (continued) 
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Figure B.9 (continued) 
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Figure B.9 (continued) 
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Figure B.10. Results for bulk temperature variation with plate spacing (eb vs. L) for CO. 
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Figure B.10 (continued) 






B-74 



& 


& 


au 


& 


Figure B.10 (continued) 
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Figure B.10 (continued) 
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Figure B.10 (continued) 
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Figure B.10 (continued) 
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Figure B.10 (continued) 
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Figure B.ll. Results for bulk temperature variation with plate spacing (eb vs. 1) for C02- 
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Figure B.ll (continued) 
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Figure B.ll (continued) 
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Figure B.ll (continued) 








Figure B.ll (continued) 
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Figure B.ll (continued) 
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Figure B.ll (continued) 
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Figure B.ll (continued) 
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Figure B.ll (continued) 
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Figure B.12 Results for bulk temperature variation with plate spacing (e b vs. L) for H 2 0 
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Figure B.12 (continued) 
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Figure B.12 (continued) 
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Figure B.12 (continued) 
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Figure B.12 (continued) 
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APPENDIX C 


COMPUTER CODE FOR TRANSIENT ENERGY TRANSFER IN 
INCOMPRESSIBLE LAMINAR FLOWS 

Following the numerical procedure discussed under the section "Method 
of Solution," a general computer program is developed to calculate the 
temperature variation and bulk temperature for a fully-developed laminar 
flow between two parallel plates having uniform wall heat flux. The 
program provides the general solution as well as the limiting (optically 
thin and large path length) solutions. The program is written for the 
three-band C0 2 gas, but can be used for any gas by employing appropriate 
thermophysical and radiative properties. 



ooo oooooooooooooooooo 


PROGRAM TRAN (0UTPUT=65,TAPE6, TAPE1 1,TAPE15, TAPE 10, TAPE1,TAPE2, 


&TAPE3) 

c*********************************************************************** 


THIS PROGRAM CALCULATES THE TRANSIENT RADIATIVE 
INTERACTION IN FULLY DEVELOPED LAMINAR FLOW OF 
C02 BETWEEN TWO PARALLEL BLACK WALLS WITH CONSTANT 
HEAT FLUX. 

IT IS BASED ON REPORT BY DR. S.N. TIWARI (REPORT NO 1-NAG 432) 


L DISTANCE BETWEEN PLATES (IN CM) 

TW WALL TEMPRATURE ( K) 

P PRESSURE (ATM) 

TT TIME (NON DIMENSIONAL) 

Z NON DIMENSIONAL DISTANCE BETWEEN THE PLATES 
Z = 0.0 AT LOWER WALL 
Z = 1.0 AT UPPER WALL 

EXTERNAL F211,F212,F213,F221,F222,F223,F511,F512,F513,F611,F612, 
&F613,F521,F522,F523,F621,F622,F623 
REAL L,KFB 

DIMENSION U (60) , EPS ( 3) , PRES (10) , TEMP (8) ,ZZ(60) ,TT (60) ,TB(65) 
&,TB0(65) ,TBL(65) ,DZZ(20) 

COMMON F1,F2,F3,B1,B2,B3,Z 

DATA U/0. 1,0. 2,0. 4, 0.5, 1.0, 2., 3., 4., 6., 7., 8., 10., 20., 40., 60., 80. 
&, 100., 200., 500., 1000./ 

DATA PRES/0.01,0.1,1.0,10.0/ 

DATA TEMP/ 500 . 0 , 1000 . 0 , 2000 . 0 , 3000 . 0/ 

DATA TT/0.05,0.1,0.5,1.0/ 

DATA DZZ/0.2,0.3,0.4, .5, .6, .7, .8/ 

WRITE (6, 110) 

WRITE(10,120) 

WRITE ( 15 , 130) 

NP=20 

NT=4 


TIME MARCHING 

DO 22 ITT=1,NT 
T=TT(ITT) 

IB=1 

C SET IB=0 FOR CALCULTING BULK TEMPRATURE, OTHERWISE IB=0 
IM=1 

C SET IM=1 FOR QUARTIC AND IM=2 FOR QUADRATIC 
IF(IB.EQ.l) GOTO 108 
WRITE(6, 150) 

WRITE(10,150) 

WRITE(15,150) 

GOTO 109 

108 WRITE(6,118) 

WRITE(10,118) 

109 WRITE(15,118) 
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WRITE(6,140)T 
WRITE(15,140) T 
WRITE(10,140) T 
00 44 IT=1,3 
TW=TEMP( IT) 

00 55 KK=3,3 
P=PRES(KK) 

DO 66 1=1, NP 
L=U( I) 

N 1 1=7 
DZ=0. 1 

DO 77 IZ=1, Nil 
Z=DZZ( IZ) 

CALCULATION OF PLANCK FUNCTION AND ITS DERIVATIVE 

WNB BAND CENTER (1/CM) 

HCK CONSTANT (K CM) 

CCC C1*C2 (ERG-K-CM**3/SEC) 

PFDBI PLANCK FUNCTION DERIVATIVE FOR I BAND 

CESS=TW**2 
STU=TW**0.5 
SNT=TW/273.0 
HCK=1. 439257246 
CCC=0. 000053847734 

SPECTROSCOPIC PROPERTIES OF C02 

HERE WE HAVE CONSIDERED ONLY THREE BANDS(15,4.3 & 2.7 MICRONS) 
WNB I BAND CENTER (1/CM) 

WNB1=667.0 
WNB2=2347.0 
WNB3=3716.0 
C2B1=HCK*WNB1 
C2B2=HCK*WNB2 
C2B3=HCK*WNB3 
CCC=0. 000053847734 
CCB1=CCC*(WNB1**4) 

CCB2=CCC*(WNB2**4) 

CCB3=CCC*(WNB3**4) 

TB1=C2B1/TW 
TB2=C2B2/TW 
TB3=C2B3/TW 
TEB1=EXP(TB1) 

TEB2=EXP(TB2) 

TEB3=EXP(TB3) 

SNTB1=CESS*( (TEB 1-1.0) **2.0) 

SNTB2=CESS*((TEB2-1.0)**2.0) 

SNTB3=CESS*((TEB3-1.0)**2.0) 

PFDB1=(CCB1*TEB1)/SNTB1 
PFDB2=(CCB2*TEB2)/SNTB2 
PFDB3=(CCB3*TEB3)/SNTB3 
C 

C BAND MOOEL CORRELATIONS (TIEN & LOWDER WIDE BAND MODEL) 
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C AZI AO I (1/CM) 

C CZSI CO 1**2 (1/ (ATM-CM) ) 

C BSI B**2 (NON DIMENSIONAL) 

C OMGI WAVE NUMBER( 1/CM) 

C SI INTEGRATED BAND INTENSITY (1/(ATM CM**2)) 

C 

AKl=(TW/300.0)**0. 5 
AK2= ( 300 . O/TW ) ** 1 . 5 
AZ1=1.29*STU 
AZ2=1.15*STU 
AZ3=2.4*STU 
C 

C KFB THERMAL CONDUCTIVITY (ERG/CM-SEC-K) 

C 

KFB=( 1488. 36517 1)*(SNT**1. 23) 

0MG1=1351.0 

0MG2=667.0 

0MG3=2396.0 

TX=-(HCK/TW) 

TX1=TX*0MG1 

TX2=TX*0MG2 

TX3=TX*0MG3 

ETX1=EXP(TX1) 

ETX2=EXP(TX2) 

ETX3=EXP(TX3) 

Cl=l . 0-ETX1 
C2=1.0-ETX2 
C3=1.0-ETX3 
BRKT=TX*( OMG 1+OMG 3 ) 

PHI2=(1.0-EXP(BRKT) )/(Cl*C3) 

TS=300.0/TW 
S1=339.685*TS 
S2=2702. 7*TS 
S3=71.497*TS*PHI2 
CZS1=S1/AZ1 
CZS2=S2/AZ2 
CZS3=S3/AZ3 
C 

C PL PRESSURE PATH LENGTH (ATM-CM) 

C 

PL=P*L 

UZ1=CZS1*PL 

UZ2=CZS2*PL 

UZ3=CZS3*PL 

BS1=0.0841*AK1 

8S2=0.32895*AK1 

PHI3=1.0+0.053*((TW/100.0)**1.5) 

DE L2= ( PH 1 3**2 - 0 ) / (PHI2*AK1) 

BS3=0.1112*DEL2 

C 

C PEI EFFECTIVE PRESSURE FOR EACH BAND (NON DIMENSIONAL) 
C 

PE1=(1.3*P)**0.7 

PE2=(1.3*P)**0.8 

PE3=(1.3*P)**0.65 
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8ETAI LINE STRUCTURE PARAMETER 

BETA1=BS1*PE1 
BETA2=BS2*PE2 
BETA3=BS3*PE3 

CORRELATION FOR EACH BAND 

F 1=2. 94* ( 1.0-EXP (- (2. 60*BETA1 ) ) ) 

F2=2.94*( 1.0-EXP (-( 2. 60*BETA2) )) 
F3=2.94*(1.0-EXP(-(2.60*BETA3))) 

B1=1.5*UZ1 
B2=1.5*UZ2 
B3=1.5*UZ3 
H1=AZ1*PFDB1 
H2=AZ2*PFDB2 
H3=AZ3*PFDB3 

12 = 1*1 

CZ=1./(Z-Z2) 

NUMERICAL INTEGRATION 

FOR DETAILS OF INTEGRATION REFER MATHEMATICAL LIBRARY 
AT NASA LANGLEY (DOCUMENT N-3) 

EPS(1)=1E-12 
EPS(2)=1E-12 
X=0.0 

IF(IM.EQ.l) GOTO 678 
CALL CADRE (X,Z,F2 11, EPS, ITEXT ,G211,IERR) 

CALL CADRE (X,Z,F212, EPS, ITEXT ,G212,IERR) 

CALL CADRE (X,Z,F2 13, EPS, ITEXT ,G213,IERR) 

CALL CADRE (Z, 1. ,F221 .EPS, ITEXT ,G221,IERR) 

CALL CADRE ( Z , 1 . , F222 , EPS , ITEXT ,G222 , IERR) 

CALL CADRE(Z,1.,F223, EPS, ITEXT.G223, IERR) 

SUM21= (G21 1-G221) 

SUM22=(G212-G222) 

SUM23=(G213-G223) 

BR1=H1*UZ1*SUM21 
BR2=H2*UZ2*SUM22 
BR3=H3*UZ3*SUM23 
BR=BR1+BR2+BR3 
AJ2=CZ* ( 2+1 . 5*L/KFB*BR ) 

GT2=12* ( EXP ( -AJ2*T) - 1 ) /AJ2 
THETA=GT2*(Z-Z2) 

GOTO 679 
678 CONTINUE 

CALL CADRE(X,Z,F511, EPS, ITEXT.G511, IERR) 

CALL CADRE(X, Z.F512.EPS, ITEXT, G512, IERR) 

CALL CADRE(X,Z,F513, EPS, ITEXT.G513, IERR) 

CALL CADRE(Z,1., F521.EPS, ITEXT, G521, IERR) 

CALL CADRE ( Z , 1 . , F522 , EPS , ITEXT ,G522 , IERR) 

CALL CADRE(Z,1. , F523.EPS, ITEXT, G523, IERR) 

CALL CADRE(X,Z, F611.EPS, ITEXT, G611, IERR) 

CALL CADRE(X, Z,F612, EPS, ITEXT.G612, IERR) 
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CALL CADRE (X,Z, F613.EPS, ITEXT ,G613,IERR) 

CALL CADRE ( Z , 1 . , F621 , EPS , ITEXT ,G621 , IERR) 

CALL CADRE(Z,1.,F622, EPS, ITEXT.G622, IERR) 

CALL CADRE (Z , 1 . , F623 , EPS , ITEXT ,G623 , I ERR) 

SUM51=G51 1-G521 

SUM52=G512-G522 

SUM53=G513-G523 

SUM61=G611-G621 

SUM62=G612-G622 

SUM63=G613-G623 

BR51=H1*UZ1*SUM51 

BR52=H2*UZ2*SUM52 

BR53=H3*UZ3*SUM53 

BR61=H1*UZ1*SUM61 

BR62=H2*UZ2*SUM62 

BR63=H3*UZ3*SUM63 

XZ=Z-2*Z**3.+Z**4. 

YZ=Z**2-Z 

ZZ=12*(Z-Z**2) 

AJ5=ZZ+1 . 5*L/KFB*(BR51+BR52+BR53) 
AJ6=-2.+1.5*L/KFB*(BR61+BR62+BR63) 

DAL=AJ5*T/XZ 

DAL1=AJ6*T/YZ 

DS-EXP(-DAL) 

DS1=EXP(-DAL1) 

GT=ZZ/(2.*AJ5)*(DS-1) 

HT=ZZ/(2.*AJ6)*(DS1-1) 

THETA=GT*XZ+HT*YZ 
679 TB(IZ)=THETA*(Z-Z2) 

IF(IB.EQ.l) GOTO 207 
WRITE(6,160) TW,P,L, THETA, Z 
207 CONTINUE 

LARGE PATH LENGTH 

IF(IM.EQ.l) GOTO 801 
DS=ABS(Z/(Z-1) ) 

AJ=AL0G(DS)+AL0G(-2/Z+l/Z**2+l)*Z+2 
H=H1+H2+H3 

AJ L=CZ* ( 2+H*L*A J/KFB ) 

GL= 12* ( EXP ( -A J L*T ) - 1 ) /AJ L 
THEL=GL*(Z-Z2) 

GOTO 802 
801 CONTINUE 
Z1=Z-1 

ALl=Zl**3./3.-3.*Z*Zl**2./2.+3.*Z*Z*Zl-Z**3.*AL0G(l-Z)-Z**3. 
&*(ll./6.-AL0G(Z) ) 

AL2=-(Zl**2./2.-2.*Z*Zl+Z**2*AL0G(l-Z)+Z*Z*(1.5-AL0G(Z) )) 
AL3=-(AL0G(1-Z)-AL0G(Z) ) 

AL=4.*AL1-6.*AL2+AL3 

AJ5=ZZ+L/KFB*AL*(H1+H2+H3) 

ALL1=-1.-Z*AL0G(1.-Z)+Z*AL0G(Z) 

ALL2=AL3 

ALL=2.*ALL1-ALL2 

AJ6=-2+L/KFB*ALL*(Hl+H2+H3) 
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DAL=AJ5*T/XZ 

DAL1=AJ6*T/YZ 

DS=EXP(-DAL) 

DS1=EXP(-DAL1) 

GT=ZZ/(2.*AJ5)*(DS-1) 

HT=ZZ/(2.*AJ6)*(DS1-1) 

THEL=GT*XZ+HT*YZ 
802 TBL( IZ)=THEL*(Z-Z2) 

IF(IB.EQ.l) GOTO 208 
WRITE(10,160) TW,P,L,THEL,Z 

208 CONTINUE 
C 

C OPTICALLY THIN 
C 

IF(IM.EQ.l) GOTO 701 
AN1=H1*UZ1 
AN2=H2*UZ2 
AN3=H3*UZ3 

AN=L/KFB*(AN1+AN2+AN3) 

AJ0=2*CZ+3*AN 

G0= 12+ ( EXP ( -AJO+T) - 1 ) /AJO 

THE0=G0*(Z-Z2) 

GOTO 702 

701 CONTINUE 

AN1=P*L*L/KFB*(S1*PFDB1+S2*PFDB2+S3+PFDB3) 

EE=Z**4-2+Z**3+Z 

FF=Z**2-Z 

AJ5=ZZ+3.*AN1*EE 

AJ6=-2.+3.*ANl*FF 

DAL=AJ5*T/XZ 

DAL1=AJ6*T/YZ 

DS=EXP(-DAL) 

DS1=EXP(-DAL1) 

GT=ZZ/(2.*AJ5)*(DS-1) 

HT=ZZ/(2.*AJ6)*(0S1-1) 

THEO=GT*XZ+HT*YZ 

702 TB0(IZ)=THE0*(Z-Z2) 

IF(IB.EQ.l) GOTO 209 
TB0(IZ)=THE0*(Z-Z2) 

IF(IB.EQ.l) GOTO 209 
WRITE(15,160) TW,P,L,THEO,Z 

209 CONTINUE 

77 CONTINUE 

IF(IB.EQ.O) GOTO 234 
TBUL=0. 

TBULL=0. 

TBUL0=0. 

DO 10 JB-l.Nll 
TBUL=TBUL+TB(JB) 

TBULL=TBULL+TBL(JB) 

10 TBU L0=TBU LO+TBO ( J B ) 

BULK=6.*DZ*TBUL 
BULKL=6.*DZ*TBULL 
BULK0=6.*DZ+TBUL0 
WRITE(6,160) TW,P,L,BULK 



WRITE ( 10 , 160) TW,P,L,BULKL 
WRITE ( 15 , 160) TW,P,L,BULKO 
234 CONTINUE 

66 CONTINUE 

55 CONTINUE 

44 CONTINUE 

33 CONTINUE 

22 CONTINUE 

110 F0RMAT(7X,*GENERAL SOLUTION FOR C02 (CONSTANT FLUX CASE)*) 

1 18 FORMAT ( 4X , *TW* , T18 ,*P* ,T28 , *L* , T38 ,*TBULK*) 

120 FORMAT (7X,*LARGE PATH LENGTH FOR C02 (CONSTANT FLUX CASE)*) 

130 F0RMAT(7X,*0PTICALLY THIN FOR C02 (CONSTANT FLUX CASE)*) 
140 F0RMAT(10X,*TIME=*,F8.3) 

150 F0RMAT(4X,*TW*,T18,*P*,T28,*L*,T38,*THETA*,T48,*Z*) 

160 F0RMAT(7(1X, F9.4)) 

STOP 

END 

FUNCTION F21 1 (ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B1*(Z-ZB) 

DENl=(Fl*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*Fl) 

AUD1=(F1*( (XX**2.0)+(4.0*XX*F1)+(4.0*F1)))/DEN1 

F211=(1-2*ZB)*AUD1 

RETURN 

ENO 

FUNCTION F212(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B2*(Z-ZB) 

DEN2=(F2*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F2) 

AU02=(F2*((XX**2.0)+(4.0*XX*F2)+(4.0*F2)))/DEN2 

F212=(1-2*ZB)*AUD2 

RETURN 

ENO 

FUNCTION F213(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B3*(Z-ZB) 

DEN3=(F3*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F3) 

AU03*(F3*( (XX**2,0)+(4.0*XX*F3)+(4.0*F3)))/DEN3 

F213=(1-2*ZB)*AUD3 

RETURN 

END 

FUNCTION F221 ( ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B1*(ZB-Z) 

DEN1=(F1*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F1) 

AUD1=(F1*( (XX**2.0)+(4.0*XX*F1)+(4.0*F1)))/DEN1 

F221=(1-2*ZB)*AUD1 

RETURN 

ENO 

FUNCTION F222(ZB) 

COMMON Fl t F2,F3,Bl,B2,B3,Z 
XX=B2*(ZB-Z) 

DEN2=(F2*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F2) 

AUD2=(F2*((XX**2.0)+(4.0*XX*F2)+(4.0*F2)))/DEN2 

F222=(1-2*ZB)*AUD2 



RETURN 

END 

FUNCTION F223(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B3*(ZB-Z) 

DEN3=(F3*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F3) 

AU03=(F3*( (XX**2.0)+(4.0*XX*F3)+(4.0*F3)))/DEN3 

F223=(1-2*ZB)*AUD3 

RETURN 

END 

FUNCTION F51 1 ( ZB ) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B1*(Z-ZB) 

DEN1=(F1*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*Fl) 

AUD1=(F1*( (XX**2.0)+(4.0*XX*F1)+(4.0*F1)))/DEN1 

F511=( 1-6*ZB*ZB+4*ZB**3)*AUD1 

RETURN 

END 

FUNCTION F512 ( ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B2*(Z-ZB) 

DEN2=(F2*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F2) 

AUD2=(F2*( (XX**2.0)+(4.0*XX*F2)+(4.0*F2)))/DEN2 

F512=(1-6*ZB*ZB+4*ZB**3)*AUD2 

RETURN 

END 

FUNCTION F513(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B3*(Z-ZB) 

DEN3=(F3*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F3) 

AUD3=(F3*((XX**2.0)+(4.0*XX*F3)+(4.0*F3)))/DEN3 

F513=(1-6*ZB*ZB+4*ZB**3)*AUD3 

RETURN 

END 

FUNCTION F521 (ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B1*(ZB-Z) 

DEN1=(F1*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*Fl) 

AUD1=(F1*((XX**2.0)+(4.0*XX*F1)+(4.0*F1)))/DEN1 

F521=(1-6*ZB*ZB+4*ZB**3)*AUD1 

RETURN 

END 

FUNCTION F522(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B2*(ZB-Z) 

DEN2=(F2*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F2) 

AUD2=(F2*( (XX**2.0)+(4.0*XX*F2)+(4.0*F2)) )/DEN2 

F522= ( 1-6*ZB*ZB+4*ZB**3)*AUD2 

RETURN 

END 

FUNCTION F523(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B3*(ZB-Z) 

DEN3=(F3*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F3) 
AUD3=(F3*( (XX**2.0)+(4.0*XX*F3)+(4.0*F3)))/DEN3 
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F523=(1-6*ZB*ZB+4*ZB**3)*AUD3 

RETURN 

END 

FUNCTION F611(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B1*(Z-ZB) 

DEN1=(F1*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*Fl) 

AUD1=(F1*( (XX**2.0)+(4.0*XX*F1)+(4.0*F1) ) )/DENl 

F61 1= ( 2*ZB- 1 ) *AUD 1 

RETURN 

END 

FUNCTION F612(ZB) 

COMMON F1,F2,F3,B1 ,B2,B3,Z 
XX=B2*(Z-ZB) 

DEN2=(F2*( (XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F2) 

AUD2=(F2*( (XX**2.0)+(4.0*XX*F2)+(4.0*F2) ) )/DEN2 

F612=(2*ZB-1)*AUD2 

RETURN 

END 

FUNCTION F613(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B3*(Z-ZB) 

DEN3=(F3*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F3) 

AUD3=(F3*( (XX**2.0)+(4.0*XX*F3)+(4.0*F3)))/DEN3 

F613=(2*ZB-1)*AUD3 

RETURN 

END 

FUNCTION F621(ZB) 

COMMON FI,F2,F3,B1,B2,B3,Z 
XX=B1*(ZB-Z) 

DENl=(Fl*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*Fl) 

AUD1=(F1*((XX**2.0)+(4.0*XX*F1)+(4.0*F1)))/DEN1 

F621=(2*ZB-I)*AUD1 

RETURN 

END 

FUNCTION F622(ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B2*(ZB-Z) 

DEN2=(F2*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F2) 

AUD2=(F2*((XX**2.0)+(4.0*XX*F2)+(4.0*F2)))/DEN2 

F622=(2*ZB-1)*AUD2 

RETURN 

END 

FUNCTION F623 ( ZB) 

COMMON F1,F2,F3,B1,B2,B3,Z 
XX=B3*(ZB-Z) 

DEN3=(F3*((XX**2.0)+(2.0*XX)+2.0)+XX)*(XX+2.0*F3) 

AUD3=(F3*((XX**2.0)+(4.0*XX*F3)+(4.0*F3)))/DEN3 

F623=(2*ZB-1)*AUD3 

RETURN 

END 



